
*For Margaret dropbox folders 
global NCREIF "/Users/BeckaBrolinson/Dropbox/NCREIF/data" 
global build 	"$NCREIF/build" 
global analysis "$NCREIF/analysis"
global results 	"$analysis/results"
global figures 	"$analysis/figures" 

*------------------------------------------------------------------------------*
*	Step 5- Event Study of Rent and utility Expenditure - unmatched	   *
*------------------------------------------------------------------------------*
{
	*Event Study 
	set more off 
	use "$build/01_annualizeddatawcontrol.dta", replace
	
	gen sqft1000s = sqft / 1000 
	gen HDD1000s = HDD / 1000 
	gen CDD1000s = CDD / 1000 
	gen used_space1000s = used_space/ 1000
	gen real_capex_ti_sqft = real_capex_ti / sqft 
	gen real_capex_bldimp_sqft = real_capex_bldimp / sqft 
	
	*List of covariates 
	#d ;
	global Cov0 	""; 
	global Cov1 	"Covered_E real_elecprice real_gasprice Unemployment HDD1000s CDD1000s 
	age age2 used_space1000s
	 percentleased real_capex_ti_sqft real_capex_bldimp_sqft "; 
	global Cov2 	"Covered_E real_elecprice real_gasprice Unemployment HDD1000s CDD1000s
	 age age2 percentleased used_space1000s
	dFundType1 dFundType2 dFundType4 dFundType5 dFundType6";
	global Cov3	"Covered_E real_elecprice Unemployment HDD1000s CDD1000s 
	age age2 used_space1000s"; 
	#d cr 
	
	*Generate Variable for Intensity of Treatment 
	bysort propnum: gen rate=everrated if year>=firstyearrated
	*Generate different bins for rating 	
	gen rate1=1 if rate>=75 & rate<82 
	replace rate1=0 if missing(rate1)
	gen rate2=1 if rate>=82 & rate<91 
	replace rate2=0 if missing(rate2)
	gen rate3=1 if rate>=91 & rate<100 
	replace rate3=0 if missing(rate3)
	*Generate a catogory for these so that I can plot them 
	gen rate_cat=. 
	replace rate_cat=1 if rate1==1 
	replace rate_cat=2 if rate2==1 
	replace rate_cat=3 if rate3==1 
	
	label var interaction "Cert*Post" 
	label var treat "1[treat=1]"
	label var logrealrentpersf "ln(Rent/Sq. Ft.) (\\$)"
	label var Covered_E "Benchmarking Law" 
	label var real_elecprice "Avg. Elec. Price (\\$/MWh)" 
	label var real_gasprice "Avg. Gas Price (\\$/Mcf)" 
	label var Unemployment "Unemployment" 
	label var age "Building Age" 
	label var age2 "Building Age Squared" 
	label var used_space "Used Space (Pct. Leased * Sq. Ft.)"
	label var used_space1000s "Used Space (Pct. Leased * Sq. Ft.) (1000s)"
	labe var percentleased "Percent Leased (\%)"
	label var real_capex_ti "Cap Exp. (Tenant Imp.)" 
	label var real_capex_ti_sqft "Cap Exp./Sq. Ft. (Tenant Imp.) (\\$)"
	label var real_capex_bldimp "Cap Exp. (Bldg. Imp.)"
	label var real_capex_bldimp_sqft "Cap Exp./Sq. Ft. (Bldg. Imp.) (\\$)"
	label var HDD "HDD" 
	label var HDD1000s "HDD (1000s)" 
	label var CDD "CDD" 
	label var CDD1000s "CDD (1000s)" 
	label var rl_yr_rentpersf "Rent per Sq. Ft."
	label var logrealutilpersf "ln(Util/Sq. Ft.) (\\$)"
	label var rl_yr_rentpersf "Rent per Sq. Ft."
	label var rl_yr_utilpersf "Util. per Sq. Ft."
	label var yrbuilt "Year Built" 
	label var sqft "Square Feet"
	label var sqft1000s "Square Feet (1000s)" 
	
	
	*Generate time since Treatment 
	gen years_pre_post_cert= year-firstyearrated
	*Generate indicators for 15 years pre and post cert (30 indicators)
	tostring years_pre_post_cert, replace
	replace years_pre_post_cert=subinstr(years_pre_post_cert, "-", "_", 1 )
	local years "_15 _14 _13 _12 _11 _10 _9 _8 _7 _6 _5 _4 _3 _2 _1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15"
	foreach year of local years{
	generate pre`year'=1 if years_pre_post_cert=="`year'"
	}
	
	*Generate the interaction variable: 
	local years "_15 _14 _13 _12 _11 _10 _9 _8 _7 _6 _5 _4 _3 _2 _1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15"
	foreach year of local years{
	generate interaction`year'=treat*pre`year'
	replace interaction`year'=0 if missing(interaction`year')
	label variable interaction`year' "`year'"
	}
	*Generate the interaction variable with the intensity var : 
	local years "0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15"
	local intensity "1 2 3" 
	foreach year of local years{
	foreach i of local intensity{
	generate interaction`i'_`year'=rate`i'*pre`year'
	replace interaction`i'_`year'=0 if missing(interaction`i'_`year')
	label variable interaction`i'_`year' "Rating Bin `i' `year'"
	}
	}	

	#d ;
	global interaction "interaction_15 interaction_14 interaction_13 interaction_12
	interaction_11 interaction_10 interaction_9 interaction_8
	interaction_7 interaction_6 interaction_5 interaction_4
	interaction_3 interaction_1 interaction0
	interaction1 interaction2 interaction3 interaction4 
	interaction5 interaction6 interaction7 interaction8 
	interaction9 interaction10 interaction11 interaction12
	interaction13 interaction14 interaction15";
	global rate1_interaction "interaction_15 interaction_14 interaction_13 interaction_12
	interaction_11 interaction_10 interaction_9 interaction_8
	interaction_7 interaction_6 interaction_5 interaction_4
	 interaction_3 interaction_1 interaction1_0
	interaction1_1 interaction1_2 interaction1_3 interaction1_4 
	interaction1_5 interaction1_6 interaction1_7 interaction1_8 
	interaction1_9 interaction1_10 interaction1_11 interaction1_12
	interaction1_13 interaction1_14 interaction1_15";
	
	global rate2_interaction "interaction2_0 interaction2_1 interaction2_2 interaction2_3 interaction2_4 
	interaction2_5 interaction2_6 interaction2_7 interaction2_8 
	interaction2_9 interaction2_10 interaction2_11 interaction2_12
	interaction2_13 interaction2_14 interaction2_15";
	
	global rate3_interaction "interaction3_0 interaction3_1 interaction3_2 interaction3_3 interaction3_4 
	interaction3_5 interaction3_6 interaction3_7 interaction3_8 
	interaction3_9 interaction3_10 interaction3_11 interaction3_12
	interaction3_13 interaction3_14 interaction3_15";
	#d cr 
	
	*Fe regression of util Exp with Heteroscedasticity and autocorrelation-consistent asymptotic variance (HAC)
	preserve
	eststo UnmatchedUtilEventStudy: parmby "reghdfe logrealutilpersf $Cov1 rl_yr_rentpersf $interaction, absorb(i.propnum i.city_cat#i.year) vce(cluster cbsa) ", label norestore
	*findit eclplot 
	*Plot the point estimates 
	keep if parmseq>=19 & parmseq<=33
	*Change the _ to - in the label variable 
	replace label=subinstr(label, "_", "-", 1 )
	sencode label if parm!="_cons", gene(parmlab)

	*Plot the outcomes 
	eclplot estimate min95 max95 parmlab, ytitle("Coefficient Estimates") ///
	xtitle("Treatment Effect By Years Before/After Certification") yline(0, lcolor(navy*.5)) ///
	xlabel(#19) bgcolor(white) graphregion(color(white)) ylab(,nogrid) legend(on lab(1 "95% CI")) 
	gr export "$figures/08_util_eventstudy_unmatched.png", replace	
	restore
	
	*Fe regression with Heteroscedasticity and autocorrelation-consistent asymptotic variance (HAC)
	preserve
	eststo UnmatchedRentEventStudy: parmby "reghdfe logrealrentpersf $Cov1 $interaction, absorb(i.propnum i.city_cat#i.year) vce(cluster cbsa) ", label norestore
	*findit eclplot 
	*Plot the point estimates 
	keep if parmseq>=18 & parmseq<=37
	*Change the _ to - in the label variable 
	replace label=subinstr(label, "_", "-", 1 )
	sencode label if parm!="_cons", gene(parmlab)

	*Plot the outcomes 
	eclplot estimate min95 max95 parmlab, ytitle("Coefficient Estimates") ///
	xtitle("Treatment Effect By Years Before/After Certification") yline(0, lcolor(navy*.5)) ///
	xlabel(#19) bgcolor(white) graphregion(color(white)) ylab(,nogrid) legend(on lab(1 "95% CI")) 
	gr export "$figures/09_util_eventstudy_unmatched.png", replace
	restore	

	*re-run to put estimates into latex table 
	eststo UnmatchedUtilEventStudy: reghdfe logrealutilpersf $Cov1 rl_yr_rentpersf $interaction, absorb(i.propnum i.city_cat#i.year) vce(cluster cbsa)
	estfe  UnmatchedUtilEventStudy, labels(year "Year FE" city_cat#year "City by Year FE" propnum "Property FE")
	eststo UnmatchedRentEventStudy: reghdfe logrealrentpersf $Cov1 $interaction, absorb(i.propnum i.city_cat#i.year) vce(cluster cbsa)
	estfe  UnmatchedRentEventStudy, labels(year "Year FE" city_cat#year "City by Year FE" propnum "Property FE")

		
	label var interaction_5 "5 Years Pre-Cert"
	label var interaction_4 "4 Years Pre-Cert"
	label var interaction_3 "3 Years Pre-Cert"
	label var interaction_1 "1 Years Pre-Cert"
	label var interaction0 "Certification Year"
	label var interaction1 "1 Year Post-Cert"
	label var interaction2 "2 Year Post-Cert"
	label var interaction3 "3 Year Post-Cert"
	label var interaction4 "4 Year Post-Cert"
	label var interaction5 "5 Year Post-Cert"
	
	
	*make Latex Table 
	esttab UnmatchedRentEventStudy UnmatchedUtilEventStudy using "$results/03_UnmatchedEventStudy.tex" , label replace booktabs ///
	alignment(SS) ///
	b(%12.3f) se(%12.3f) star(* 0.05) /// sets format of parameters, standard errors, and stars
	stats(N, fmt(%9.0fc)) /// adds comma to Observation number
	nonotes nogaps /// removes notes from bottom of table 
	indicate(`r(indicate_fe)', label("Y" "")) /// adds y & N for FE inclusions
	drop(interaction_15 interaction_14 interaction_13 interaction_12 /// 
	interaction_11 interaction_10 interaction_9 interaction_8 ///
	interaction_7 interaction_6 interaction6 interaction7 interaction8 ///
	interaction9 interaction10 interaction11 interaction12 ///
	interaction13 interaction14 interaction15 ) /// dropping the first and last 10 interaction terms to be able to fit table on one page
	title(Unmatched Sample: Energy Star Certification Event Study \label{tabappendix01eventstudyunmatched}) ///
	addnote({\scriptsize * p$<$0.05. The standard errors reported in parenthesis have been clustered at the CBSA level.})
	
	
	*Event Study 
	set more off 
	use "$build/01_annualizeddatawcontrol.dta", replace
	*Generate Variable for Intensity of Treatment 
	bysort propnum: gen rate=everrated if year>=firstyearrated
	*Generate different bins for rating 	
	gen rate1=1 if rate>=75 & rate<82 
	replace rate1=0 if missing(rate1)
	gen rate2=1 if rate>=82 & rate<91 
	replace rate2=0 if missing(rate2)
	gen rate3=1 if rate>=91 & rate<100 
	replace rate3=0 if missing(rate3)
	*Generate a catogory for these so that I can plot them 
	gen rate_cat=. 
	replace rate_cat=1 if rate1==1 
	replace rate_cat=2 if rate2==1 
	replace rate_cat=3 if rate3==1 
	
	
	gen sqft1000s = sqft / 1000 
	gen HDD1000s = HDD / 1000 
	gen CDD1000s = CDD / 1000 
	gen used_space1000s = used_space/ 1000
	gen real_capex_ti_sqft = real_capex_ti / sqft 
	gen real_capex_bldimp_sqft = real_capex_bldimp / sqft 
	
	set more off 
	*List of covariates 
	#d ;
	global Cov0 	""; 
	global Cov1 	"Covered_E real_elecprice real_gasprice Unemployment HDD1000s CDD1000s 
	age age2 used_space1000s
	 percentleased real_capex_ti_sqft real_capex_bldimp_sqft "; 
	global Cov2 	"Covered_E real_elecprice real_gasprice Unemployment HDD1000s CDD1000s
	 age age2 percentleased used_space1000s
	dFundType1 dFundType2 dFundType4 dFundType5 dFundType6";
	global Cov3	"Covered_E real_elecprice Unemployment HDD1000s CDD1000s 
	age age2 used_space1000s"; 
	#d cr 

	label var interaction "Cert*Post" 
	label var treat "1[treat=1]"
	label var logrealrentpersf "ln(Rent/Sq. Ft.) (\\$)"
	label var Covered_E "Benchmarking Law" 
	label var real_elecprice "Avg. Elec. Price (\\$/MWh)" 
	label var real_gasprice "Avg. Gas Price (\\$/Mcf)" 
	label var Unemployment "Unemployment" 
	label var age "Building Age" 
	label var age2 "Building Age Squared" 
	label var used_space "Used Space (Pct. Leased * Sq. Ft.)"
	label var used_space1000s "Used Space (Pct. Leased * Sq. Ft.) (1000s)"
	labe var percentleased "Percent Leased (\%)"
	label var real_capex_ti "Cap Exp. (Tenant Imp.)" 
	label var real_capex_ti_sqft "Cap Exp./Sq. Ft. (Tenant Imp.) (\\$)"
	label var real_capex_bldimp "Cap Exp. (Bldg. Imp.)"
	label var real_capex_bldimp_sqft "Cap Exp./Sq. Ft. (Bldg. Imp.) (\\$)"
	label var HDD "HDD" 
	label var HDD1000s "HDD (1000s)" 
	label var CDD "CDD" 
	label var CDD1000s "CDD (1000s)" 
	label var rl_yr_rentpersf "Rent per Sq. Ft."
	label var logrealutilpersf "ln(Util/Sq. Ft.) (\\$)"
	label var rl_yr_rentpersf "Rent per Sq. Ft."
	label var rl_yr_utilpersf "Util. per Sq. Ft."
	label var yrbuilt "Year Built" 
	label var sqft "Square Feet"
	label var sqft1000s "Square Feet (1000s)" 
	
	
	
	*Generate time since Treatment 
	gen years_pre_post_cert= year-firstyearrated
	*Generate indicators for 15 years pre and post cert (30 indicators)
	tostring years_pre_post_cert, replace
	replace years_pre_post_cert=subinstr(years_pre_post_cert, "-", "_", 1 )
	local years "_15 _14 _13 _12 _11 _10 _9 _8 _7 _6 _5 _4 _3 _2 _1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15"
	foreach year of local years{
	generate pre`year'=1 if years_pre_post_cert=="`year'"
	}
	
	*Generate the interaction variable: 
	local years "_15 _14 _13 _12 _11 _10 _9 _8 _7 _6 _5 _4 _3 _2 _1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15"
	foreach year of local years{
	generate interaction`year'=treat*pre`year'
	replace interaction`year'=0 if missing(interaction`year')
	label variable interaction`year' "`year'"
	}
	*Generate the interaction variable with the intensity var : 
	local years "0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15"
	local intensity "1 2 3" 
	foreach year of local years{
	foreach i of local intensity{
	generate interaction`i'_`year'=rate`i'*pre`year'
	replace interaction`i'_`year'=0 if missing(interaction`i'_`year')
	label variable interaction`i'_`year' "Rating Bin `i' `year'"
	}
	}	

	#d ;
	global interaction "interaction_15 interaction_14 interaction_13 interaction_12
	interaction_11 interaction_10 interaction_9 interaction_8
	interaction_7 interaction_6 interaction_5 interaction_4
	interaction_3 interaction_1 interaction0
	interaction1 interaction2 interaction3 interaction4 
	interaction5 interaction6 interaction7 interaction8 
	interaction9 interaction10 interaction11 interaction12
	interaction13 interaction14 interaction15";
	global rate1_interaction "interaction_15 interaction_14 interaction_13 interaction_12
	interaction_11 interaction_10 interaction_9 interaction_8
	interaction_7 interaction_6 interaction_5 interaction_4
	 interaction_3 interaction_1 interaction1_0
	interaction1_1 interaction1_2 interaction1_3 interaction1_4 
	interaction1_5 interaction1_6 interaction1_7 interaction1_8 
	interaction1_9 interaction1_10 interaction1_11 interaction1_12
	interaction1_13 interaction1_14 interaction1_15";
	
	global rate2_interaction "interaction2_0 interaction2_1 interaction2_2 interaction2_3 interaction2_4 
	interaction2_5 interaction2_6 interaction2_7 interaction2_8 
	interaction2_9 interaction2_10 interaction2_11 interaction2_12
	interaction2_13 interaction2_14 interaction2_15";
	
	global rate3_interaction "interaction3_0 interaction3_1 interaction3_2 interaction3_3 interaction3_4 
	interaction3_5 interaction3_6 interaction3_7 interaction3_8 
	interaction3_9 interaction3_10 interaction3_11 interaction3_12
	interaction3_13 interaction3_14 interaction3_15";
	#d cr 
	
	
	*Same estimation as above, but including 3 different treatment intensities 
	preserve
	parmby "reghdfe logrealutilpersf $Cov1  $rate1_interaction $rate2_interaction $rate3_interaction, absorb(i.propnum i.city_cat#i.year) vce(cluster cbsa) ", label norestore
	
	replace label=subinstr(label, "_", "-", 1 )
	replace label=subinstr(label, "Rating Bin 1", "", .) 
	replace label=subinstr(label, "Rating Bin 2", "", .)
	replace label=subinstr(label, "Rating Bin 3", "", .)
	sencode label if parm!="_cons", gene(parmlab)
	
	*Keep only TE parameter Estimates 
	keep if parmseq>=18
	drop if parmseq>=38 & parmseq<=42 | parmseq>=54 & parmseq<=58 | parmseq>=70 & parmseq<=74 
	*relabel some data to prepare to plot 
	gen time=label
	destring time, replace 
	
	twoway 	(connected estimate time if parmseq>=17 & parmseq<=25, lcolor(navy)) ///
	(line min95 time if parmseq>=17 & parmseq<=25, lcolor(gs10) lpattern(dot)) /// 
	(line max95 time if parmseq>=17 & parmseq<=25, lcolor(gs10) lpattern(dot)) ///
	(connected estimate time if parmseq>=26 & parmseq<=35, lcolor(dkgreen*.5) mcolor(dkgreen*.5)) /// 
	(line min95 time if parmseq>=26 & parmseq<=35, lcolor(gs10) lpattern(dash)) /// 
	(line max95 time if parmseq>=26 & parmseq<=35, lcolor(gs10) lpattern(dash)) ///	
	(connected estimate time if parmseq>=42 & parmseq<=52, lcolor(dkgreen*.75) mcolor(dkgreen*.75)) /// 
	(line min95 time if parmseq>=42 & parmseq<=52, lcolor(gs10) lpattern(shortdash)) /// 
	(line max95 time if parmseq>=42 & parmseq<=52, lcolor(gs10) lpattern(shortdash)  yline(0, lcolor(navy*.5)))	///	
	(connected estimate time if parmseq>=52, lcolor(dkgreen) mcolor(dkgreen)) /// 
	(line min95 time if parmseq>=52, lcolor(gs4) lpattern(dot)) /// 
	(line max95 time if parmseq>=52, lcolor(gs4) lpattern(dot) ///
	bgcolor(white) graphregion(color(white)) ///
	title("Intensity of Certification Effect") xlabel(#19) ///
	legend(order(1 "Pre-Treatment" 2 "Pre-Treatment CI" 4 "Rating (75-82)" ///
	5 "Rating1 CI" 7 "Rating (82-90)" 8 "Rating 2 CI" 10 "Rating (91-100)" ///
	11 "Rating 3 CI"  )))	
	gr export "$figures/10_util_eventstudy_unmatched_intensity3.png", replace	

	twoway 	(connected estimate time if parmseq>=17 & parmseq<=25, lcolor(navy)) ///
	(line min95 time if parmseq>=17 & parmseq<=25, lcolor(gs10) lpattern(dot)) /// 
	(line max95 time if parmseq>=17 & parmseq<=25, lcolor(gs10) lpattern(dot)) ///
	(connected estimate time if parmseq>=26 & parmseq<=35, lcolor(dkgreen*.5) mcolor(dkgreen*.5)) /// 
	(line min95 time if parmseq>=26 & parmseq<=35, lcolor(gs10) lpattern(dash)) /// 
	(line max95 time if parmseq>=26 & parmseq<=35, lcolor(gs10) lpattern(dash)) ///	
	(connected estimate time if parmseq>=42 & parmseq<=52, lcolor(dkgreen*.75) mcolor(dkgreen*.75)) /// 
	(line min95 time if parmseq>=42 & parmseq<=52, lcolor(gs10) lpattern(shortdash)) /// 
	(line max95 time if parmseq>=42 & parmseq<=52, lcolor(gs10) lpattern(shortdash)  yline(0, lcolor(navy*.5))	///	
	bgcolor(white) graphregion(color(white)) ///
	title("Intensity of Certification Effect") xlabel(#19) ///
	legend(order(1 "Pre-Treatment" 2 "Pre-Treatment CI" 4 "Rating (75-82)" ///
	5 "Rating1 CI" 7 "Rating (82-90)" 8 "Rating 2 CI" )))	
	gr export "$figures/11_util_eventstudy_unmatched_intensity2.png", replace	
	restore 

*Rather than looking at the above (by time treated) plot, how about TE by year? 
	set more off 
	use "$build/01_annualizeddatawcontrol.dta", replace
		
	gen sqft1000s = sqft / 1000 
	gen HDD1000s = HDD / 1000 
	gen CDD1000s = CDD / 1000 
	gen used_space1000s = used_space/ 1000
	gen real_capex_ti_sqft = real_capex_ti / sqft 
	gen real_capex_bldimp_sqft = real_capex_bldimp / sqft 
	
	
	#d ; 
	local Years "2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012
	2013 2014 2015"; 
	#d cr 
	foreach year of local Years{
	*can only assess for first year treated because by including the outcome variable 
	*in the years before treatment means that a building cannot have been treated 
	*in years before or else unconfoundedness doesn't hold 
	gen rate`year'=1 if year==firstyearrated & firstyearrated==`year' 
	replace rate`year'=0 if rate`year'!=1 
	}
	
	
*Generate the interaction variable: 
	#d ; 
	local Years "2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012
	2013 2014 2015"; 
	#d cr 
	foreach year of local Years{
	generate year`year'cert=treat*rate`year'
	replace year`year'cert=0 if missing(year`year'cert)
	label variable year`year'cert "`year'"
	}
#d ; 
	global annualcert  "year2001cert year2002cert year2003cert 
	year2004cert year2005cert year2006cert year2007cert 
	year2008cert year2009cert year2010cert year2011cert 
	year2012cert year2013cert year2014cert year2015cert"; 
#d cr 

*Estimate Annual Effects 
preserve
	*we are readjusting years since treatment for every time a building gets recertified 
	*Relative to year 2000 
	parmby "reghdfe logrealutilpersf $Cov1 rl_yr_rentpersf $annualcert,  absorb(i.propnum i.city_cat#i.year) vce(cluster city_cat)", label norestore
	*findit eclplot 
	*Plot the point estimates 
	keep if parmseq>=14 & parmseq<=28
	*Change the _ to - in the label variable 
	sencode label if parm!="_cons", gene(parmlab)

	*Plot the outcomes 
	eclplot estimate min95 max95 parmlab, ytitle("Coefficient Estimates") ///
	xtitle("Years Before and After First Certification") yline(0, lcolor(navy*.5)) ///
	xlabel(#19, angle(30)) bgcolor(white) graphregion(color(white)) ylab(,nogrid) legend(on lab(1 "95% CI")) 
	gr export "$figures/12_util_eventstudy_unmatched_year.png", replace
	restore
	
	*Same as above but for rent 
	preserve 
	parmby "reghdfe logrealrentpersf $Cov1 $annualcert,  absorb(i.propnum i.city_cat#i.year) vce(cluster city_cat)", label norestore
	keep if parmseq>=13 & parmseq<=27
	*Change the _ to - in the label variable 
	sencode label if parm!="_cons", gene(parmlab)

	*Plot the outcomes 
	eclplot estimate min95 max95 parmlab, ytitle("Coefficient Estimates") ///
	xtitle("Years Before and After First Certification") yline(0, lcolor(navy*.5)) ///
	xlabel(#19, angle(30)) bgcolor(white) graphregion(color(white)) ylab(,nogrid) legend(on lab(1 "95% CI")) 
	gr export "$figures/13_rent_eventstudy_unmatched_year.png", replace
}


*Heterogenous effects by energy star score 
{
*load the unmatched sample to do event study with het effects 
	
	set more off 
	use "$build/01_annualizeddatawcontrol.dta", replace
		
	gen sqft1000s = sqft / 1000 
	gen HDD1000s = HDD / 1000 
	gen CDD1000s = CDD / 1000 
	gen used_space1000s = used_space/ 1000
	gen real_capex_ti_sqft = real_capex_ti / sqft 
	gen real_capex_bldimp_sqft = real_capex_bldimp / sqft 
	
	
	*Generate Variable for Intensity of Treatment 
	bysort propnum: gen rate=everrated if year>=firstyearrated
	*Generate different bins for rating 	
	gen rate1=1 if rate>=75 & rate<81 
	replace rate1=0 if missing(rate1)
	gen rate2=1 if rate>=81 & rate<100 
	replace rate2=0 if missing(rate2)
	
	*Generate a catogory for these so that I can plot them 
	gen rate_cat=. 
	replace rate_cat=1 if rate1==1 
	replace rate_cat=2 if rate2==1 
	
*Diff in diff regression with Het treatment effects 

	set more off 
	*List of covariates 
	#d ;
	global Cov0 	""; 
	global Cov1 	"Covered_E real_elecprice real_gasprice Unemployment HDD1000s CDD1000s 
	age age2 used_space1000s
	 percentleased real_capex_ti_sqft real_capex_bldimp_sqft "; 
	global Cov2 	"Covered_E real_elecprice real_gasprice Unemployment HDD1000s CDD1000s
	 age age2 percentleased used_space1000s
	dFundType1 dFundType2 dFundType4 dFundType5 dFundType6";
	global Cov3	"Covered_E real_elecprice Unemployment HDD1000s CDD1000s 
	age age2 used_space1000s"; 
	#d cr 
	

	*Generate time since Treatment 
	gen years_pre_post_cert= year-firstyearrated
	*Generate indicators for 15 years pre and post cert (30 indicators)
	tostring years_pre_post_cert, replace
	replace years_pre_post_cert=subinstr(years_pre_post_cert, "-", "_", 1 )
	local years "_15 _14 _13 _12 _11 _10 _9 _8 _7 _6 _5 _4 _3 _2 _1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15"
	foreach year of local years{
	generate pre`year'=1 if years_pre_post_cert=="`year'"
	}
	
	*Generate the interaction variable: 
	local years "_15 _14 _13 _12 _11 _10 _9 _8 _7 _6 _5 _4 _3 _2 _1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15"
	foreach year of local years{
	generate interaction`year'=treat*pre`year'
	replace interaction`year'=0 if missing(interaction`year')
	label variable interaction`year' "`year'"
	}
	*Generate the interaction variable with the intensity var : 
	local years "0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15"
	local intensity "1 2" 
	foreach year of local years{
	foreach i of local intensity{
	generate interaction`i'_`year'=rate`i'*pre`year'
	replace interaction`i'_`year'=0 if missing(interaction`i'_`year')
	label variable interaction`i'_`year' "Rating Bin `i' `year'"
	}
	}	
	
	#d ;
	global interaction "interaction_15 interaction_14 interaction_13 interaction_12
	interaction_11 interaction_10 interaction_9 interaction_8
	interaction_7 interaction_6 interaction_5 interaction_4
	interaction_3 interaction_1 interaction0
	interaction1 interaction2 interaction3 interaction4 
	interaction5 interaction6 interaction7 interaction8 
	interaction9 interaction10 interaction11 interaction12
	interaction13 interaction14 interaction15";
	
	global rate1_interaction "interaction_15 interaction_14 interaction_13 interaction_12
	interaction_11 interaction_10 interaction_9 interaction_8
	interaction_7 interaction_6 interaction_5 interaction_4
	 interaction_3 interaction_1 interaction1_0
	interaction1_1 interaction1_2 interaction1_3 interaction1_4 
	interaction1_5 interaction1_6 interaction1_7 interaction1_8 
	interaction1_9 interaction1_10 interaction1_11 interaction1_12
	interaction1_13 interaction1_14 interaction1_15";
	
	global rate2_interaction "interaction2_0 interaction2_1 interaction2_2 interaction2_3 interaction2_4 
	interaction2_5 interaction2_6 interaction2_7 interaction2_8 
	interaction2_9 interaction2_10 interaction2_11 interaction2_12
	interaction2_13 interaction2_14 interaction2_15";
	

	#d cr 

	*differences in differences regression with two different treatment levels 
	
	*Event Study with two different treatment leves 
	
	preserve
	parmby "reghdfe logrealutilpersf $Cov1 rl_yr_rentpersf $rate1_interaction $rate2_interaction , absorb(i.propnum i.city_cat#i.year) vce(cluster cbsa) ", label norestore
	
	replace label=subinstr(label, "_", "-", 1 )
	replace label=subinstr(label, "Rating Bin 1", "", .) 
	replace label=subinstr(label, "Rating Bin 2", "", .)

	
	sencode label if parm!="_cons", gene(parmlab)

	*Keep only TE parameter Estimates 
	keep if parmseq>=14
	drop if parmseq > 68
	*relabel some data to prepare to plot 
	gen time=label
	destring time, replace 
	

	twoway 	(connected estimate time if parmseq>=19 & parmseq<=27, lcolor(navy)) ///
	(line min95 time if parmseq>=19 & parmseq<=27, lcolor(gs10) lpattern(dot)) /// 
	(line max95 time if parmseq>=19 & parmseq<=27, lcolor(gs10) lpattern(dot)) ///
	(connected estimate time if parmseq>=28 & parmseq<=38, lcolor(dkgreen*.5) mcolor(dkgreen*.5)) /// 
	(line min95 time if parmseq>=28 & parmseq<=38, lcolor(gs10) lpattern(dash)) /// 
	(line max95 time if parmseq>=28 & parmseq<=38, lcolor(gs10) lpattern(dash)) ///	
	(connected estimate time if parmseq>=44 & parmseq<=54, lcolor(dkgreen*.75) mcolor(dkgreen*.75)) /// 
	(line min95 time if parmseq>=44 & parmseq<=54, lcolor(gs10) lpattern(shortdash)) /// 
	(line max95 time if parmseq>=44 & parmseq<=54, lcolor(gs10) lpattern(shortdash)  yline(0, lcolor(navy*.5))	///	
	bgcolor(white) graphregion(color(white)) ///
	xlabel(#19) ylab(,nogrid) ///
	legend(order(1 "Pre-Treatment" 2 "Pre-Treatment CI" 4 "Rating (75-80)" ///
	5 "Rating1 CI" 7 "Rating (81-100)" 8 "Rating 2 CI" )))	
	gr export "$figures/15_utilhetestar_eventstudy_unmatch.png", replace	
	restore 
	
	preserve
	parmby "reghdfe logrealrentpersf $Cov1 $rate1_interaction $rate2_interaction , absorb(i.propnum i.city_cat#i.year) vce(cluster cbsa) ", label norestore
	
	replace label=subinstr(label, "_", "-", 1 )
	replace label=subinstr(label, "Rating Bin 1", "", .) 
	replace label=subinstr(label, "Rating Bin 2", "", .)

	
	sencode label if parm!="_cons", gene(parmlab)

	*Keep only TE parameter Estimates 
	keep if parmseq>=18
	drop if parmseq > 68
	*relabel some data to prepare to plot 
	gen time=label
	destring time, replace 
	

	twoway 	(connected estimate time if parmseq>=18 & parmseq<=26, lcolor(navy)) ///
	(line min95 time if parmseq>=18 & parmseq<=26, lcolor(gs10) lpattern(dot)) /// 
	(line max95 time if parmseq>=18 & parmseq<=26, lcolor(gs10) lpattern(dot)) ///
	(connected estimate time if parmseq>=27 & parmseq<=37, lcolor(dkgreen*.5) mcolor(dkgreen*.5)) /// 
	(line min95 time if parmseq>=27 & parmseq<=37, lcolor(gs10) lpattern(dash)) /// 
	(line max95 time if parmseq>=27 & parmseq<=37, lcolor(gs10) lpattern(dash)) ///	
	(connected estimate time if parmseq>=43 & parmseq<=53, lcolor(dkgreen*.75) mcolor(dkgreen*.75)) /// 
	(line min95 time if parmseq>=43 & parmseq<=53, lcolor(gs10) lpattern(shortdash)) /// 
	(line max95 time if parmseq>=43 & parmseq<=53, lcolor(gs10) lpattern(shortdash)  yline(0, lcolor(navy*.5))	///	
	bgcolor(white) graphregion(color(white)) ///
	xlabel(#19) ylab(,nogrid) ///
	legend(order(1 "Pre-Treatment" 2 "Pre-Treatment CI" 4 "Rating (75-80)" ///
	5 "Rating 1 CI" 7 "Rating (81-100)" 8 "Rating 2 CI" )))	
	gr export "$figures/15_renthetestar_eventstudy_unmatch.png", replace	
	restore 
	
	}
